#load gplot2 library
library(ggplot2)
#Read data 
library(readxl)
Data1 <- read_excel("Data1.xlsx", sheet = "Leaf width_day 9 & 30")
View(Data1)

#convert categorical variables to factors 
str(Data1)
Data1$EMS = as.factor(Data1$EMS)
Data1$Accession = as.factor(Data1$Accession)
Data1$DAP = as.factor(Data1$DAP)
Data1$Leaf_width = as.numeric (Data1$Leaf_width)
str(Data1)

#two-way ANOVA with no interaction 
two.way <- aov(Leaf_width ~ EMS + Accession + DAP, data = Data1)
#two-way ANOVA with interaction
interaction <- aov(Leaf_width ~ EMS * Accession * DAP, data = Data1)

#AIC calculates the best-fit model by finding the model that explains the 
#largest amount of variation in the response variable while using the fewest 
#parameters (Bevans, 2023). 
#https://www.scribbr.com/statistics/two-way-anova/#:~:text=A%20two%2Dway%20ANOVA%20is,combination%2C%20affect%20a%20dependent%20variable.
#install the AIC package
install.packages("AICcmodavg")
#Load the AIC package
library(AICcmodavg)
#Run the AIC package
model.set <- list(two.way, interaction)
model.names <- c("two.way", "interaction")

aictab(model.set, modnames = model.names)

#Summary of two-way ANOVA
summary(two.way)

#Check for homoscedasticity. 
#Check whether the model fits the assumption of homoscedasticity
par(mfrow=c(2,2))
plot(two.way)
par(mfrow=c(1,1))

#Post-hoc test. 
#Tukey’s Honestly-Significant-Difference (TukeyHSD) test
turkey.two.way<-TukeyHSD(two.way)
turkey.two.way
#plot of turkey test
tukey.plot.aov<-aov(Leaf_width ~ EMS:Accession:DAP, data = Data1)
tukey.plot.test<-TukeyHSD(tukey.plot.aov)
tukey.plot.test
plot(tukey.plot.test, las = 1)

#Presentation of result in bar chart
Hypocotyl_length_plot <- ggplot(Data1, aes(x = Accession, 
                                           y = Hypocotyl_length,
                                           fill= EMS, 
                                           
                                           
) 
) +
  
  
  stat_summary(fun = "mean", geom = "bar", 
               position = "dodge", 
               colour = "gray25", 
               alpha = .7)+
  stat_summary(fun = "mean", 
               geom = "point",
               position = position_dodge(0.9),
               size = 1)+
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  
  stat_summary(fun.data = 'mean_se', geom = 'errorbar', na.rm = TRUE,  
               position = position_dodge(0.9), width = 0.4, color = 'black')+
  theme(
    strip.text.x = element_text(margin = margin(2, 0, 2, 0))
  )+
  ylim(0, 7) +
  scale_fill_grey()
Hypocotyl_length_plot
Hypocotyl_length_plot <- Hypocotyl_length_plot +
  labs(title = "Effect of EMS treatment on hypocotyl length of Bambara groundnut",
       x = "Bambara groundnut accessions",
       y = "Hypocotyl length (cm)") +
  theme(axis.title = element_text(color = "black",
                                  face='bold',
                                  size = 12), 
        legend.title=element_text(color = "black", 
                                  face = "bold", 
                                  size = 12),
        legend.text = element_text(color = "black", 
                                   face = "bold", 
                                   size = 12), 
        plot.title = element_text(color = "black", 
                                  face = "bold", 
                                  size = 12), 
        axis.text = element_text(color = "black",
                                 face = "bold",
                                 size = 8))
Hypocotyl_length_plot
ggsave("Hypocotyl length_DAP.png")